% PAIRANALYSIS is a collection of codes for the analysis of paired
% recordings data. This is the starting point for the analysis. Based on
% the following codes, I am going to write more specific functions.

% written by Taro Kiritani, tarokiritani2008@u.northwestern.edu
% 10/15/20010

rect = getrect;

% check the pair is connected or not.
% 3sd of noise is the threshold




% plot ppr and amplitude
data = [get(get(gco,'Children'),'XData');get(get(gco,'Children'),'YData')];

response_begin = find(data(1,:)>0);
response_end = find(data(1,:)<0.2);

baseline = intersect(response_begin,response_end);
baseline = data(2,baseline);

epsp.first.amp = max(data(2,[2000:2400]) - mean(baseline));
epsp.second.amp = max(data(2,[2500:2900]) - mean(baseline));
epsp.third.amp = max(data(2,[3000:3400]) - mean(baseline));
epsp.fourth.amp = max(data(2,[3500:3900]) - mean(baseline));


ppr = epsp.second.amp/epsp.first.amp
figure(h)
plot(epsp.first.amp,ppr,'Marker','*','Color','r');hold on

figure(h2)
plot([1,epsp.second.amp/epsp.first.amp,epsp.third.amp/epsp.first.amp,epsp.fourth.amp/epsp.first.amp],'Color','r');hold on


% quantal analysis.
directoryName = 'C:\DATA\Taro\CELLS\TK0350';
files = dir(directoryName);
str = {files.name};
[s,v] = listdlg('PromptString', 'Select files:', 'SelectionMode','multiple','ListString',str);


for i = 1:length(str(s))
    fullfilename{i} = fullfile(directoryName, str{s(i)});
    [pathstr{i}, tracename, ext, versn] = fileparts(fullfilename{i});
    data.completename = fullfile(pathstr{i}, [tracename, ext, versn]);
    % open the file, extract the data
    currentTrace = load(data.completename, '-mat');
%     currentTrace.userFcn.filename{i} = data.completename{i};
    baseline = mean(currentTrace.data.ephys.trace_2(1:2000));
    epsp(i) = max(currentTrace.data.ephys.trace_2(2000:2500))
end

figure;hist(epsp-baseline)

% extract baseline-subtracted data.
h = findobj(get(gca, 'Children'), 'Type', 'line')
d = cell2mat(get(h, 'YData'));
baseline = mean(d(:,1:1999),2)
barray = repmat(baseline, 1, 20000);
bsarray = d - barray;
figure;plot(bsarray')

% normalization with the fits epsp.
firstepsp = bsarray(:,2000:2500);
firstepsp = max(firstepsp,[],2);
figure;plot((bsarray./repmat(firstepsp,1,20000))');

% epsp train amplitude analysis. measure the epsp from the trough. intended
% to look at fascilitation/depression.
base{1} = dnorm(:,2450:2500);
base{2} = dnorm(:,2950:3000);
base{3} = dnorm(:,3450:3500);
% base{4} is for the recovery pulse.
base{4} = dnorm(:,8450:8500);

for k = 1:3
% base{k} = mean(base{k},2);
epsp(:,k) = max(dnorm(:,(2500 + 500*(k-1)):(3000 + 500*(k-1))),[],2);
end
% base{4} = mean(base{4},2);
% epsp(:,4) = min(dnorm(:,8500:9000),[],2) - base{4};
epsp(:,4) = max(dnorm(:,8500:8600),[],2) - mean(dnorm(:,8450:8500),2);

% epsp train amplitude analysis. measure the trough from the baseline. intended
% to look at integration.
base = dnorm(:,1950:2000);
base = mean(base,2);
for k = 1:3
    trough{k} = dnorm(:,2450 + 500 * (k-1):2500 + 500 * (k-1));
    trough{k} = mean(trough{k},2);
%     trough{k} = trough{k} - base;
end
figure;plot(cell2mat(trough)')



% analysis for the decay of epsps
% select axes and run the next scripts.
linehundles = get(gco,'Children');
figure
for m = 1:length(linehundles)
data = [get(linehundles(m),'XData');get(linehundles(m),'YData')];
xdata = intersect(find(data(1,:) > 0.87),find(data(1,:) < 0.94));
ydata = data(2,xdata);
% baseline subtraction
baseline = data(2,(intersect(find(data(1,:) > 0.94),find(data(1,:) < 0.95))));
ydata = ydata - mean(baseline);

xdata = data(1,xdata);
[estimates, model] = fitExp(xdata-xdata(1),ydata);
f = @(t) estimates(1) * exp(-estimates(2) * (t - xdata(1)));
plot(xdata,ydata);hold on
ezplot(f, [0.87 0.94])
timeConspn(m) = 1/estimates(2);
end


% plot in subfigure(4,4,13) in epsc.fig
% get data in subfigure(445)
h = findobj(get(gca, 'Children'), 'Type', 'line');
d = cell2mat(get(h, 'YData'));
for n = 1:4
    baseline{n} = mean(d(:,1950 + 500 * (n - 1):2000 + 500 * (n - 1)),2);
    epsc{n} = min(d(:,2000 + 500*(n - 1):2500 + 500*(n - 1)), [], 2);
end
baseline{5} = mean(d(:,8450:8500),2);
epsc{5} = min(d(:,8500:9000), [], 2);
Baseline = cell2mat(baseline);
EPSC = cell2mat(epsc);
EPSC = EPSC - Baseline;